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Abstract 

In two-phase flow simulations, a difficult issue is usually the treatment 
of surface tension effects. These cause a pressure jump that is proportional 
to the curvature of the interface separating the two fluids. Since the 
evaluation of the curvature incorporates second derivatives, it is prone 
to numerical instabilities. Within this work, the interface is described 
by a level-set method based on a discontinuous Galerkin discretization. 
In order to stabilize the evaluation of the curvature, a patch-recovery 
operation is employed. There are numerous ways in which this filtering 
operation can be applied in the whole process of curvature computation. 
Therefore, an extensive numerical study is performed to identify optimal 
settings for the patch-recovery operations with respect to computational 
cost and accuracy. 


1 Introduction and motivating example 

In order to simulate immiscible two-phase flows, one usually has to consider 
surface tension effects. These induce a pressure jump which is proportional to 
the total curvature k of the fluid interface 3. Precisely, the momentum equation 
for the fluid domains 21 and fB is given as 

d 

— (pit) + div (pu (g) It) + Vt/; = pAit, (1) 

ot 
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while at the fluid interface 3 = 21 fl *8 the velocity u and pressure i/' are coupled 
via the jump condition 

|(v^/ - ^ {Vu H- Vii^)) na] = (2) 

see e.g. 0 or 0 . We briefly introduce the notation required for this formula: 

• The computational domain ft C is decomposed into fluid 21, fluid *8 
and the {D — l)-dimensional interface U, i.e. = 21 U J U 18. Regarding 
this work, we restrict ourselves to the case D = 2. 

• denotes the normal vector on 3, oriented so that “it points from 21 to 
18” and k denotes the total curvature of 3. We assume 3 to be smooth 
enough, so that both, and k are at least in C°(3). 

• the jump operator for u S C°(n \ 3) is defined as 

|u] (x) = lim {u{x + - u{x - ^nj)). (3) 


• /i and p denote the dynamic viscosity and the density of the fluid, which are 
usually constant within either 21 and 18, but have a jump at the interface. 

This setting may be described by a level-set function p € which fulfills 


< 0 

in 

21 


p> Q 

in 

<8 

( 4 ) 

p = Q 

on 

3 



Then, obviously, 

n -3 = and k = div (nj). (5) 

\Vp\2 

The latter expression is also called Bonnet’s formula. In dependence of the 
level-set gradient Vp and the level-set Hessian d'^p, it can be expressed as 


with 


ViVy’l; 


curv(g, H) = 


=: curv(Vi^, d'^p) 

( 6 ) 

tr(ff) g^Hg 

\g \2 \g\l 

( 7 ) 


Note that by the introduction of p, the properties and k, which were initially 
only defined on 3, were smoothly extended to the whole domain H. 

The purpose of this paper is to benchmark various hltering strategies for 
Bonnet’s formula, based on patch recovery. In order to assess the quality of 
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a filtering strategy, the L^-error in the curvature may only be one property of 
interest; in addition, we define the following benchmark problem: 



|'0] = UK on fJ 

jS/'ih ■ n^l = 0 on fJ 

t/j = 0 on dn 


( 8 ) 


This Poisson problem shares some behavior with a two-phase flow problem o, 
since it would correspond e.g. to the pressure-projection step of a projection 
method. It is solved by an extended discontinuous Galerkin method (XDG), 
details are given in section 

1.1 Outline for this work 

The main motivation for this work is given in section 11.41 where it is demon¬ 
strated, for a simple example, that different choices for the level-set-field ip, 
with equal zero-sets J, can lead to very different results if just a straightforward 
projection of the curvature is used. Next, the patch-recovery operation which 
we employ to overcome this issues is introduced (section [5]) and the test setup is 
defined Isection [XT]) . Results and a discussion, together with recommendations 
for the use of patch-recover are given in section [3.21 Finally, in appendix El the 
discretization of the Poisson problem (|S]) is given. 

1.2 Treatment of surface tension in numerical methods. 

In multi-phase flow simulations, surface tension plays an important role as soon 
as length scales become small enough, e.g. in the range of centimeters and below. 
Over the previous decade, a huge variety of two-phase flow solvers have been 
proposed and a significant portion of them uses level-set methods. Some of them 
use a smeared interface approach - where the piecewice constant properties p 
and p, in eqs. (ITI2]) are regularized so that their derivatives are finite - like 
e.g. the works of Herrmann Q. In contrast to this, extended finite element or 
discontinuous Galerkin methods use an approximation space that is conformal 
with the location of the discontinuities, like e.g. done by Gross and Reusken [^. 
Common to all this methods is that surface tension forces have to be computed 
from a level-set field p. 

From a flow-solver perspective, what finally matters is the force which is 
induced into the momentum balance m by surface tension effects. This force 
is only active at the interface T itself, i.e. in some finite volume (FV), finite 
element (FE) or discontinuous Galerkin (DG) discretizations, it usually appears 
at the right-hand-side of the momentum equation as 




or 


Here, v denotes a test function (in 2D FV-methods usually, v = (1, 0) and v = 
(0,1)) while 5o denotes a delta-distribution on T. It is either singular, thereby 
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both integrals given above are equivalent, or regularized, i.e. ‘smeared out’ in 
a small neighborhood around 3. Whether a ‘sharp’ or a ‘smeared’ approach 
is used usually depends on the design of the flow solver, i.e. whether the flow 
solver represents the density and viscosity jump in a ‘sharp’ or smeared fashion. 

However, the evaluation of curvature is not the only way to compute surface 
tension forces. Using a variant of Stoke’s integral theorem, one can re-express 
surface tension forces in some control volume K d as 


® (TKn 3 • u dS = ® a{I — n:i n^) : Vv dS — Y as ■ v d£, (9) 

danif Yanif JapnK) 

where s is tangential to 3 and normal onto d{3 fl K). and jf ... d^ denotes an 
integral over the boundary of the surface 3 H K; i.e. in 2D, J... d£ it is a 0- 
dimensional point integral (counting measure) over the two end-points of the 
line 3 n K, while in 3D it is the 1-dimensional line integral over the boundary 
of the surface 3 (1 K. Since Kn:i can also be expressed by the Laplace-Beltrami 

- operator on the manifold 3 - for details, we refer to Gross and Reusken 01 

- the reformulation of the surface tension is often referred to as the Laplace- 
Beltrami approximation of the surface tension terms. The main benefit of this 
re-formulation is that the right-hand-side of eq. ([9]) does not depend on second 
derivatives of ip. This seems very attractive e.g. in a FV-discretiztion, where 
the test function is constant and therefore the first term on the right-hand-side 
of ([9]) vanishes and the problem is reduced to the computation of tangents s to 
3. 

The Laplace-Beltrami formulation also shows good results in the extended 
finite element (XFEM) context, as demonstrated e.g. by Gross and Reusken 0, 
Gheng and Fries 0, Sauerland and Fries [10 . 

However, if the numerical integration over the surface 3 becomes more precise 

- as is the case in this work - equation is also fulfilled numerically and 
therefore one cannot expect significant ‘filtering properties’ from the Laplace- 
Beltrami formulation, i.e. any oscillations present in the level-set would be 
captured equally by both, the left- and the right-hand-side of equation (0. 

In the work of Chen et al. 0, a projection of the curvature on the broken 
piecewise polynomial space is obtained from the relation 





dK 


S/p 

\Vp\2 


■ riKV dS 



• Vv dV 


( 10 ) 


for a test function v and a control volume K with outer normal tik. Since this is 
equivalent to the piecewise polynomial reconstruction of Bonnet’s formula ®, 
one can expect that it is prone to the same issues as we outline in section [L4l 
However, the right-hand-side of equation [10] provides the option to use a non¬ 
local flux-formulation for at the boundary integral. This was exploited 

e.g. by Heimann 0, were an additional diffusion term is incorporated into the 
non-local projection to deal with discontinuities. 
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A completely different approach to curvature evaluation is proposed by 
Heimann et al. [6|. They use an extended DG space - they refer to it as ‘unfitted 
DG’ - for velocity and pressure, a conventional DG field for the level-set while 
for curvature computation, they employ finite difference scheme. 

The patch-recovery post-processing was first introduced for finite element 
methods by Zienkiewicz and Zhu [I 3 . To our knowledge, the only application 
of patch recovery (see section [2]) to curvature evaluation has been demonstrated 
by Herrmann and Jibben Q. There, a nodal projection is used, while this work 
employs an L^-projection, making the algorithm independent of any choice of 
nodes and more general with respect to the shape of the patches. The main 
achievement of this work is that we also investigate the application of patch- 
recovery for ‘intermediate’ values, like gradients V(p and Hessians of the 
level-set function ip. 

1.3 Piecewise polynomial approximation 

For sake of completeness, we briefly introduce the numerical grid as well as the 
approximation spaces used e.g. for (p or ip. These are fairly standard, and can 
be found in similar form in many textbooks, see e.g. Hesthaven and Warburton 
0 or Di Pietro and Ern Q. We define: 

• the numerical grid: = {ATi,..., ATj}, with h being the maximum di¬ 

ameter of all cells Kj. The cells cover the whole domain (D = IJ^ Af,), 
but do not overlap 1 dx = 0 for I ^ j). We restrict ourselves 

to non-curved grids, i.e., each cell Kj can be described as the image of 
a reference cell ATref under an affine-linear mapping Tj : —)• i.e. 

T,(Kref) = K,; 

• the broken piecewise polynomial approximation space of maximum total 
degree p: 

:= < ./ G VKGJih:flK= G 

^ 0<i+j<p 

• the broken piecewise polynomial approximation space of maximum 

p: 

Qp{fih) := >1 / G A^(D); V K £ Hh ■ /|x = x'‘y^b,j, bij G 

• the continuous piecewise polynomial space of order p: 

:= Qp(f?„)nC°(D) 
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• the extended broken approximation space of order p, which is used to 
discretize the Poisson problem ([8]): 

:=Fp{M.h)lm®Fp{Rh)l<s, (14) 

where lx denotes the characteristic function for set X C fl. 

• the projector onto some subspace X of L^(n) in the L^-norm: 

L^(f2) 9 / I—Proj^ (/) S X. (15) 

• that all derivative operators (div(—), V, A and the Hessian d^) in this 
work should be understood in a broken sense. We assert that all properties 
defined here are almost everywhere sufficiently smooth, so that the respec¬ 
tive derivatives exist in a classical sense everywhere up to a set of measure 
zero. Therefore, we do not make the usual distinctions, e.g. between the 
gradient V and the broken gradient V h ■ 

1.4 Motivating example. 

The main issue with using Bonnet’s formula is that the results are extremely 
sensitive to minor disturbances in the level-set function (^, which we are going 
to illustrate by the following example. The exact interface 3ex is chosen to be 
a circle with radius R = 0.8 around the origin. Obviously, there are infinitely 
many choices for ip. 

For a sharp-interface method, like the one described in appendix El the 
interface 3 has to be a closed surface, therefore it is required that also the 
numerical representation of the level-set field is at least continuous, i.e. p G 
C°(n). To ensure this, the exact level-set function i^ex is projected onto the 
continuous piecewise polynomial space Q 2 {^h)- 

Two representations of a circular interface around the origin with radius 
8/10 are compared, namely: 

case (a) : p = (8/10)^ — — y'^ quadratic, p = (pex 

case (b) : p = Project)(^/lO — yj signed-distance, p ^ ipex 

The domain H = (—3/2, 3/2)^ is discretized by 18 x 18 equidistant quadrilateral 
cells. Obviously, the exact curvature at the exact interface in this example 
k |3 = —10/8, and the exact solution to the Poisson problem ([5]), for a = 1/10 

_ J 0 in 21 

“ I 1/8 in 25 ■ 

Note that in case (a), the quadratic form already is a member of the space 
Q 2 (M/i), therefore the circular interface is represented exactly in the polynomial 
approximation space. In case (b), where this is not the case, the approximation 
of the circular interface still seems reasonably good: for any x in the zero-set of 
(/?, the error in the radius, i.e. | \x \2 — 8/10 | is less than 1.2 • lO”'^. 
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However, if in sub-sequence the curvature 


K := (curv(V(/ 3 , d^ip)) (17) 

is computed, and the Poisson problem ([5]) with a right-hand-side based on k is 
solved, a quite reasonable error in ip, and, more severe, is introduced, which 
can also be seen in figures [T] and [5] 


error type 

quadratic 

signed-distance 

radius 

Vx €3 : \x\2 — 8/10 

< 10-11 

< 1.2 • 10“^ 

curvature 

Va; G 3 : iK(a:)| — (—10/8) | 

< 2.2 • 10-9 

< 0.14 

pressure 

||'0 — V'exll 

< 10-1° 

Ri 1.3-10-® 

pressure gradient 

lIVV^-VV'exIL 

< 2•10"® 

« 0.4 


With respect to the two-phase Navier-Stokes equation, the error in pressure 
tp and pressure gradient V '0 are especially inadvertent: from a physical point 
of view, a circular droplet in a velocity field u = 0 represents a steady state 
with minimal energy. Since the velocity is linked to the pressure gradient by the 
momentum equation o, a non-zero pressure gradient would induce an artificial, 
non-zero velocity into a state-state solution. 

On the other hand, the quadratic test-case basically confirms that the ex¬ 
tended DG discretization of the Poisson problem is solved exactly by the nu¬ 
merics, up to machine accuracy. A necessary perquisite for this is, that that the 
quadrature rules used to implement the extended DG discretization are suffi¬ 
ciently precise, so that the discretization error is not dominated by the error of 
the numerical integration. 

So, this example demonstrates that the result of a level-set based two-phase 
computation may heavily depend on the choice of the level-set field. The purpose 
of this study is to identify filters by which this undesirable dependence can be 
minimized. 

Remark on smooth interface methods: those methods usually treat sur¬ 
face tension as a smoothed Delta-distribution, therefore they do not require 
(fi G C°(H). However, the discussion on smoothed interface methods is beyond 
the scope of this publication. 

2 Patch recovery in the — sense. 

We introduce the patch-recovery operator 

prc(^ : ^ PqiAh) ( 18 ) 

which is used to filter the level-set field ip and its derivatives. Some auxiliary 
definitions are required: 

• A cell K € Ah IS called to be cut (by the zero-level-set J) if 1 dS 7 ^ 0. 
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0.80012 


case a: quadratic 


Xi 


0.79988 



-3-2-10123 

Z{x) 


it 

-1.1 

a; 

0 

a 

> 

-1.25 

CJ 

-1.4 



-3-2-10123 

A{x) 


case b: signed-distance 


ll 

' 

M 

' 

M 

' 

y 

ll 


1 

1 

1 

1 

i 


-3-2-10123 

/:{x) 


0.80012 

0.8 

0.79988 


- 1.1 



-3-2-10123 

A{x) 


Figure 1: Polar plot of a quadratic versus a signed-distance level-set function 
for the circular interface with radius 0.8: for points a; € 9, the radius \x \2 
and the curvature at a ;, i.e. k{x) are plotted in dependence of the angular 
coordinate /(a:) For the quadratic representation, the interface is exact up to 
round-off errors, yielding a quite accurate representation of the curvature in a 
piecewise polynomial space with sufficiently high order. On the other hand, for 
the signed distance representation, an error with a magnitude of about 10“'^ in 
the radius yields an L°°-error around 10 % in the curvature. 
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0 ) 


0 ) 

> 



case (a): quadratic 
= ( 8 / 10 )^ — — y‘^ 


case (b): signed-distance 

V = (s/io - 




Figure 2: Comparison of a quadratic versus a signed-distance level-set function 
ip. While the zero level-set of both representation is pretty similar, the error in 
the curvature k for case (b) is obvious. succession, this causes high oscillations 
in the pressure gradient of the Poisson problem whose right-hand-side depends 


on K. 



• A cell L e is called to be a neighbor of cell K G if they share 
at least one point, i.e. if L f] K ^ {}. (Rem.: by defining neighbors as 
sharing one point - instead of the usual definition of sharing an edge - 
this e.g. induces that in a Cartesian grid a quadrilateral element has 8 
neighbors, and not just 4 as it would be the case with the usual definition 
of neighbors.) 

• We further define the set of all cut-cells, 

:= {K G ikh] K is cut} 
and all cut cells and their neighbors 

:= {AT G A" or one of its neighbors are cut} 
and the union of these cells, i.e. 

^cc.O _ ^cc ^cc.l _ IJ ^ 


The patch-recovery operator with width w G {0,1} is defined as the A^-projection 
onto the broken polynomial space on the composite cell Qk, which is formed 
from all neighbor cells of some cell K. Precisely: For a cell K G let 

be Li,... ,Lj G the neighbors of AT in Then one can define a 

composite cell Qk ■= A" U Ai U ... U A/ and the polynomial space of order 
q on this composite cell, IPg({(5if}). For cell K, the patch-recovery operation 
prc.j„(M) =: V is then defined as the A^-projection of u onto IPg({(5if}), i.e. 


v\k 


Pj'0jp,({QK}) (“IQk) 
0 


if a: G 

a Ki 


By prc(„(u) we denote the composition of I patch-recovery operations, and con¬ 
sequently prc° (u) = u. 


Implementation of the patch-recovery operator: For the space Pg(A/i), 
we assume an orthonormal basis {(l>j,n)j=i,...,j,n=i,...,N, with supp((()jn) = Ay. 
Here j is called the cell-index, while n is called the mode index. Then, one 
chooses a polynomial basis {0n)n=i,...,N, on the composite cell Qk and computes 
factors Anim to express the basis functions in terms of (pjn, i.e. 

^ ^ 4^limAjiim- 
i,m 

li denotes a mapping from indices i to the cell indices of the cells which make 
up Qk- To aid numerical stability, the basis may be “pre-orthonormalized”, 
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e.g. the polynomials may be chosen to be orthonormal in the bounding box of 
Qk- For the mass matrix M of one gets 


Mnm ■= / OnOm dV 
JQk 


— ^ ^ ^ ^ -^nis-^mjr j 4^lis4^ljr dV 

z,s j,r 

— ^ ^ ^nir^mir- 


—5ii 5sr 


Then, an orthonormal basis on Qk is given by 


^ ^ ^mS-nm — EE Snm-^nir j 4^11 


z,r \ n 


= .Bn 


where the change-of-basis matrix S is given by the relation 
S^MS = I, resp. M = {S-^f S-\ 


i.e. S is the inverse of the Cholesky-factor of M. Then, for some function u = 
the L^-projection of u onto is given as v := 

with 


'^m — 


U-drn dV 


Qk 

= EE*bs^^ 


j,s 




mir I 4^lir4^ljS dv 

'Qk 

V 

— Sij Sr-s 


mjs • 


Finally, v can be re-expressed in the ‘original’ basis <j)jn of the space Pq(^/i), 
within e.g. cell Ki^ as 





4>lir- 


The computationally most expensive operation is the construction of the tensors 
dlnir, Snm and finally Bmir- Once these are computed, multiple evaluations of 
the patch-recovery operator are comparatively cheap, since they can be imple¬ 
mented as matrix-vector multiplications. 
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3 The configurable curvature algorithm and Re¬ 
sults 

3.1 Test setup 

We investigate three basic cases, also shown in figure [31 these are: 

case a, large circle: := 0.8 — 

case b, small circle: (/3ex := 0.25 — \J 

case c, ’peanut’-shape: := 3 — 0.9 cos(a:) — \/{x 1)^ -I- 

-\J{x - ly + y‘^. 

For case (a) and (b), the domain is chosen as 12 = (—3/2, 3/2), discretized by 
18 X 18 equidistant cells, while for case (c) 12 = (—3, 3) x (—2, 2), discretized by 
30 X 20 equidistant cells. Since we investigate a huge variety of configurations 
for the curvature algorithm - namely 6144 configurations for each test-case, see 
below - we do not alter the DG polynomial degree of the level-set function nor 
of pressure '0 and the pressure gradient. However, the influence of polynomial 
degree of the curvature approximation n is investigated, see below. We define the 
approximation of the exact level-set field (/?ex on the broken and the continuous 
piecewise polynomial space, i.e. 

V5br := Projp^(^,^) ((^ex) (19) 

and 

ipco ■■= Pl'OjQC0(^^) (y^br). (20) 

For each test case and algorithm configuration, three different Lp' errors are 
recorded: 

curvature on cut-cells: ||k — «:ex||L 2 (Qco) 

pressure: 

pressure gradient jump: || fV-ip - Vi/exl • 'n-a||L2(a) = II [Vi/l • wa||L 2 ( 3 ) 

The latter error is especially important in the context of a two-phase Navier- 
Stokes - problem, as already noted in the first example (section 11.41) . For a 
circular interface - which represents, form a physical point of view, a natural 
minimal-energy state of a droplet - one usually wants the pressure gradient 
as low as possible. A non-zero pressure gradient would correspond to some 
‘artificial’ velocity that is introduced in the momentum equation o, preventing 
the simulation from reaching the natural, minimal-energy final state with zero 
velocity. 

Evaluation of curvature: When applying patch-recovery construction to 
curvature-evaluation, several options are at hand: Where in the algorithm 
should patch-recovery be applied? Onto the input, i.e. the level-set function, 
onto the output, i.e. the curvature, onto intermediate results like gradients or 
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case (a) case (b) case (c) 

‘large circle’ ‘small circle’ ‘peanut’ 



Figure 3: The three different test geometries which are investigated: a rather 
large circle, with respect to grid width h, a small circle and a ’peanut’-shaped 
interface. 


Hessians, or everywhere? Which polynomial degree should be chosen for the 
curvature? How many cells around the interface should be used for the patch 
recovery operation? 

In order to evaluate different configurations for patch recovery - based cur¬ 
vature algorithms, a modular algorithm with the following options was imple¬ 
mented; the filtered curvature R, used as an input to the Poisson problem (|S]), 
resp. (1^ . is computed as 


/ := one of/ 

/:= prc^(/) 

fj 

g := one of | 
g := prc^ (g) 

r d^f 

TT r 1 ^9 

H ■= one of < 

[ Vg 

H:= prc^(ff) 

K := curv ^one of ■! 
. R := prc^(K) 


9 

9 


one of 



( 21 ) 


A flowchart is given in figure ID In detail, this algorithm provides the following 
configuration options: 

• Which field should be used as an input for the algorithm (switch ‘/ =?’)? 
— the SEM-representation (case i-e. / := (pco- 
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- the DG-representation (case Vbr’)) i-®- / ■= ‘/’br- 

• From which field should the level-set gradient g be computed (switch 
‘9 =?’)? 

- from the un-filtered level-set (case ‘V/’), i.e g := V/. 

- from the filtered level-set (case ‘V/’), i.e g := V/. 

• From which field should the Hessian H of the level-set be computed 
(switch H =?’)? 

- from the un-filtered level-set (case ‘9^/’)) be. H := d'^f. 

- from the un-filtered level-set gradient (case ‘Vg’), i.e. H := Vg. 

- from the filtered level-set (case ‘9^/’), i.e. H := d^f. 

- from the filtered level-set gradient (case ‘Vg’), i.e. H := Vg. 

• The number li G {1,2,5,10} of patch-recovery-cycles that are applied to 
the level-set field and its derivatives, i.e. f prc^^(/), g := prc^^(g) and 
H :=prc'i(ff). 

• The number I 2 G {0,1, 5,10} of patch-recovery-cycles that are applied to 
the computed curvature k i.e. k := prc^^(K). 

• The DG-polynomial degree of the curvature k and the filtered derivatives 

is altered: K,K,fG ¥ap{Sih), g £ H G where p = 4 

is the polynomial degree of pdGj with the integer multiplier a G {1, 2,3}. 

3.2 Results and discussion 

For all test-cases, the three different errors measures span across several magni¬ 
tudes, see figures and [HI with one exception, namely the error W-ip — tpexWL^^n) 
in case (c), see bottom-left plot in figure [HI Therefore, the pressure error of case 
(c) is not considered in the survey. 

• It is, in all test-cases generally better to use the DG-representation (/ = 

Pbr) for curvature computation than the continuous GG (/ = pc«)- Con¬ 
sidering that the broken approximation (pbr £ is in general more 

precise than the continuous approximation G Q 2 {^h)^ this seems not 
very surprising: 
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Figure 4: Flow chart for curvature evaluation: the filtered curvature k is com¬ 
puted either from the SEM {ipc°) or the DG (:/3br) representation of the level-set 
field. The configuration of the algorithm depends on the switches for the level- 
set source (/ =?), for the gradient source {g =?), for the Hessian source [H =?), 
on the number of patch-recovery cycles for the derivatives (h) and the curvature 
itself {I 2 ), as well as the width w of the patch-recovery domain; furthermore, 
there are switches (hi. V?, fil. 9^?) to evaluate Bonnet’s formula from either 
the un-filtered gradient g or the filtered gradient g and, independently from 
either the un-filtered Hessian H or the filtered Hessian H. In addition, not 
shown in this graph, we vary the DG polynomial degree of /, g, H, k and k as 
well as the domain of the patch-recovery. 
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ly’br - ‘Pex||L2(Qcc') 

IT’CO — T’ex |L2(f2cc) 

case (a) 

3.192-10-' 

5.087-10-5 

case (b) 

1.165-10-5 

2.396 -10-4 

case (c) 

1.71 -10-5 

3.759 -10-4 


V(/3br — V(^ex (L 2 (ncc ))2 

\\V(pco - V(Pex (L2(Qcc))2 

case (a) 

3.128-10-5 

2.008 -10-5 

case (b) 

1.096-10-5 

1.096-10-2 

case (c) 

1.342 -10-5 

1.418-10-2 


— 3^(/?ex (L2(qcc'))2x2 

Wd'^ipCO - l92(^g,j ^^2(-Qcc))2x2 

case (a) 

2.173 -10-5 

9.045 -10-2 

case (b) 

7.806 -10-2 

4.364 -10-4 

case (c) 

7.958 -10-2 

4.71 - 10-4 

The minimal errors achieved for the broken and the continuous approxi 


mation, over all configurations are: 



min. over configs. with 
broken approx. 

(/ = T’br) 

min. over configs. with 
continuous approx. 

(/ = Vc«) 

curvature error ||k — b;ex||L 2 (Qcc): 

case (a) 
case (b) 
case (c) 

1.874-10-5 

3.821 -10-5 

5.46- 10-5 

8.008 -10-4 

3.606 -10-2 

6.43 - 10-2 

pressure error W'lj) - V'ex |L 2 (n) : 

case (a) 
case (b) 
case (c) 

1.053-10-5 

2.669 -10-5 
1.403-10-4 

1.697 -10-5 

6.34- 10-5 
1.542-10-4 

pressure gradient jump error 

IliVV'l-na 11^2(3) : 

case (a) 
case (b) 
case (c) 

1.928-10-4 

2.519-10-5 

5.966 -10-2 

4.19- 10-4 
3.916-10-5 
5.454-10-2 


Most interestingly, the latter error measure, || IVV:] • n 3 ||i 2 p), seems not 
to be affected by the choice of broken versus continuous approximation. 

• Regarding performance, the most influential factors are - not surprisingly 
- the polynomial order of the curvature approximation space P 4 q,(R/i), 
and the width w of the patch-recovery domain. The major computational 
cost of the patch-recovery operator is the construction of the projector 
onto the aggregate cells, while the evaluation of the patch-recovery is in 
comparison quite fast. Therefore, the number of patch-recovery cycles 
li and I 2 is only of minor influence to the run-time. For all algorithm 
configurations which employ any kind of patch-recovery, we observe the 
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following mean run-times 
different (a, w) - pairs: 


and standard deviations for 


(a, w) 

1,0 

1,1 

2,0 

2,1 

3,0 

3,1 


case (a) 


0.29 

0.7 

1.3 

6.2 

9.0 

48 

seconds 


0'tru„,c,» 

0.53 

0.84 

1.1 

2.5 

3.0 

6.9 

seconds 


f Mtrun.c,™ \ 

\ /^*run,l,0 / 

1.0 

2.5 

4.5 

22 

32 

168 

normalized 


f ^*run,(X \ 

\ ^trun.1.0 / 

1.9 

2.9 

4.0 

8.7 

10 

24 

normalized 

case (b) 


0.24 

0.38 

0.62 

2.4 

3.9 

17 

seconds 


<^trun,c,r. 

0.49 

0.62 

0.79 

1.5 

2.0 

4.1 

seconds 


( Mtrun.c,™ \ 

\ ^^run,l,0 / 

1.0 

1.6 

2.6 

9.7 

16 

71 

normalized 


( \ 

\ ^trun,l,0 / 

2.0 

2.5 

3.2 

6.3 

8.1 

17 

normalized 

case (c) 


0.25 

0.62 

1.3 

7.6 

11 

63 

seconds 



0.50 

0.79 

1.2 

2.7 

3.4 

8.0 

seconds 


\ /^*run,l,0 / 

1.0 

2.5 

5.4 

31 

46 

257 

normalized 


( '^*run,Q,iu \ 

V '"‘run, 1,0 ) 

2.0 

3.2 

4.7 

11 

14 

32 

normalized 


It becomes apparent that especially the cases with {a,w) = (3,0) and 
{a,w) = (3,1) are computationally expensive: among all test-cases, the 
{a,w) = (3,0) - configurations are between 16 and 38 times more ex¬ 
pensive than the {a,w) = (1,0) - configuration; for the {a,w) = (3,1) - 
configurations that factor is 71 and 257. 

• For all test-cases, no significant improvement int the errors Wtp — tpexW(n) 
and II(VV' — V^ex) • ^ia||L 2 ( 3 ) could be observed for a > 2 and w > 0. 
Furthermore, no significant improvement could be obtained by any kind 
of filtering of Level-Set gradient or Hessian, nor by filtering the curvature 
itself. Indeed, from figures [5] and [ 6 l it becomes obvious that configurations 
(marked as • in cited figures) which 

— use the broken polynomial-level-set field (/ = iphr), 

— perform patch recover only on cut cells itself (w = 0 ), 

— use polynomial degree 2 • p for / and k (a = 2 ), 

— compute gradient and Hessian from the filtered level-set field {g = 
V/, H = d^f), 

— employ no additional filtering of the gradient nor the Hessian (fil. V = 
false, fil. d'^ = false) and 

— employ no additional filtering of computed curvature (I 2 = 0 ) 

are among the best-performing and fastest configurations for test cases, 
with respect to each error measure. 
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• As obvious from figures [S] and [51 there are certain configurations which 
perform better in certain test cases. However, there are no configurations 
which perform significantly better overall, i.e. with respect to all three 
error measures in all three test cases. In order to check that, we selected 
the 10 fastest configurations out of the 100 best-performing, with respect 
to the three different error measures in all test-cases. From those config¬ 
urations, which performed well overall, none significantly out-performed 
the configuration recommendation given above. 

• Finally, we perform a separate discussion of those configurations which use 

the continuous representation of the level-set field (/ = While 

itself is continuous, its first and second derivatives contain higher levels of 
oscillations than the broken polynomial representation yibr of the level-set 
field, as already noted above. Here, filtering of the gradient g and the 
Hessian H provide a significant improvement. The top-5 configurations 
overall are: 


9=? 

H =? 

fil. V? 

fil. 52 ? 

a 

w 

h 

h 

V/ 

Vg 

true 

true 

1 

0 

1 

0 

V/ 

Vg 

true 

true 

1 

0 

1 

1 

V/ 

Vg 

true 

false 

1 

0 

1 

1 

V/ 

Vg 

false 

false 

1 

0 

1 

1 

V/ 

Vg 

false 

false 

1 

0 

2 

1 


Common for all of these cases is that the gradient g is compute from the 
filtered level /, the Hessian H is computed from the filtered gradient g, the 
best results were quite surprisingly achieved for lower order approximation 
spaces (i.e. 0=1, i.e. f,g,g,H,K,k are of DG-polynomial order 1 ■ a) 
and that the width of the patch-recovery domain w = 0. The performance 
of these top-5 configuration is visualized in the scatter plots figure 0 and 
111 

4 Conclusions and outlook. 

A very important result of this survey is first, that we found that there is no 
need to increase the polynomial degree of the filtered properties by more than 
a factor of two times, in comparison to the polynomial degree of the original 
level-set field; second, that we found it sufficient to perform the patch recovery 
just on the layer of cut-cells themselves, and not considering any values outside 
of that band. Since the run-time is dominated by those two factors, it is very 
pleasant that these factors can be kept low. 

Furthermore, we were able to show the additional benefit of filtering also the 
gradient and the Hessian in the case of the spectral-element representation of 
the level set. 

Within this study we focused on the curvature computation in a very con¬ 
trolled setup, where an analytic expression for the level-set field was given. 
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11 ^ ^ex||L2(r2ccj 



Legend: ^ 


configurations with / = t/Jbr, w = 0, a = 2, g = V/, 
H = d^f, fil. V = false, fil. = false, h = 0. 
best-performing configurations with respect to curvature 
error in case (a) 


Figure 5: A scatter plot, showing all algorithm configurations, for all test 
cases, presenting the runt-time of the curvature evaluation versus the error 
measure ||K-«:ex||L 2 (o cc^ . Configurations marked with symbol • are overall well¬ 
performing, i.e. they are among the most accurate and fastest configurations 
for each of the three error measures and all test cases, see also figure [H There 
are certainly configurations, like the ones marked by ♦ , which perform better 
in a specific test-case, but not overall. 
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case 


0) 

u 


0 ) 

iX! 

CJ 


0.01 

Legend: 


11^ - '*/’ex||L2(n) 



0.01 0.1 


10 


0.01 0.1 


10 


0.1 


1 


10 


100 0.01 



100 0.01 



100 0.01 


II |VV^] •na||L2(a) 



0.1 


10 


0.1 


10 


0.1 


1 


10 


100 



100 



1 

0.1 

0.01 

0.001 

0.0001 

) 

10 

1 

0.1 

0.01 

0.001 

) 

10 

1 

0.1 

0.01 


100 


configurations with / = tpbr, w = 0, a = 2, g = V/, 
H = d^f, fil. V = false, fil. = false, h = 0. 


Figure 6: A scatter plot, showing all algorithm configurations, for all test cases, 
presenting the run-time of the curvature evaluation versus the two different error 
measures Wtp — 4’ex\\L^{y) and || |Vp] ||L 2 (a). Configurations marked with symbol 
• are overall well-performing, i.e. they are among the most accurate and fastest 
configurations for each of the three error measures and all test cases, see also 
figure m 
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11 ^ Kex||z,2(f2co) 



Legend: 


• top-5 configurations overall; common configuration settings 
are: g — V/, H = Vg, a = 1, w = 0. 


Figure 7: A scatter plot, showing all algorithm configurations which use 

the continuous representation of the level-set (/ = ipco), for all test cases, 
presenting the runt-time of the curvature evaluation versus the error mea¬ 
sure \\k - «:ex||L 2 (o cc ^ . Configurations marked with symbol • are overall well¬ 
performing, i.e. they are among the most accurate and fastest configurations 
for each of the three error measures and all test cases, see also figure [51 There 
are certainly configurations, like the ones marked by x, which perform better 
in a specific test-case, but not overall. 


21 







case 


CJ 


CD 


CD 


IIV’ - V’ex||L2(o) 





0.01 0.1 1 10 100 


II IVV’l -nallL^p) 

1 

0.1 
0.01 
0.001 
0.0001 

0.01 0.1 1 10 100 

10 
1 

0.1 
0.01 
0.001 

0.01 0.1 1 10 100 

10 
1 

0.1 
0.01 

0.01 0.1 1 10 100 





Legend: 




top-5 configurations overall; common configuration settings 
features are: g — V/, H = Vg, a = 1, w = 0. 


Figure 8: A scatter plot, showing all algorithm configurations which use the 
continuous representation of the level-set (/ = for all test cases, presenting 
the run-time of the curvature evaluation versus the two different error measures 
IIV' “ V'ex||L 2 (a) and || |Vp] Hisp). Configurations marked with symbol • are 
overall well-performing, i.e. they are among the most accurate and fastest con¬ 
figurations for each of the three error measures and all test cases, see also figure 
0 
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We are currently developing a two-phase Navier-Stokes solver based on the ex¬ 
tended DG discretization presented in section El where curvature computation 
is an important part. 

This survey focused solely on patch recovery for filtering DG-based prop¬ 
erties. Obviously, there other techniques for constructing filters, e.g. WENO, 
which we may investigate in future works. 


A Extended DG discretization of the Poisson 
problem 

Finally, we give a brief description of the extended discontinuous Galerkin 
method which is used to solve the Poisson eq. (|8]). We search for {4>,v) G 
F 2 i^h) X IPs'SO that for all {f,w) G x Ff{^hY 

b{'tjj,w) = s{w) 

b{f,w) = 0 (22) 


The bilinear form &(—,—) is given as 

Kf, w) = - [ / div (w) dx- (f |t!] • {{p}} dS. 

Jn Jru3 


while the linear form s(—), which represents the surface-tension terms, is given 
as 


s{w) = j) aknj ■ |u;] dS. 

Again, V denotes also the broken gradient, and div (—) the broken divergence. 
The set T is the union of all edges of all cells, i.e. T := 

denotes an almost-everywhere continuous normal field, which is equal to on 
3, an outer normal on dVl and normal onto the cell boundaries. Given that, the 
jump- and mean-value operator are defined as 


H (x) 
{{«}}(*) 
M (x) 


lini 4 \,o (u(x + C^a.r) - u(x - Cna,r)) 
lim^N^o I (u(x + ^nr) + u(x - ^nr)) 
limjx^o u{x - ^n^,r) 
limjx^o u{x - ^na,r) 


on (r \ dV,) U fJ, 
on (r \ dfl) U 3, 
on dn, 
on dfl. 


In order to perform numerical integration on the cut cells, i.e. on domains like 
ATnSt, ATnlB, Kr\3, etc., a quadrature technique presented in [ll| is used. The 
linear system obtained from discretization (l22l) is solved by the direct sparse 

“0a 


solver PARDISO, see 15 
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